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ABSTRACT 



An algorithm for the identification of the input distribution 
matrix of a linear, stationary system operating in a stochastic 
environment is derived. The identification is accomplished by defining 
a set of autocorrelation functions for a noise element composed of a 
linear combination of the input distribution matrix elements and the 
random excitations of the system. Another possible identification 
method employing a Kalman filter is discussed and the problems associ- 
ated with its derivation are presented. Results of computer simulations 



for both methods are included. 
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I. INTRODUCTION 



In recent years, an emphasis on plant identification for a system 
operating in a stochastic environment has developed in the field of 
automatic control. It is necessary to know the dynamics of the system 
before it may be effectively controlled. In order to know the system 
dynamics, a suitable mathematical model of the system must be found. 

For a system subjected to a random forcing function, finding a suitable 
mathematical model involves the identification of both the state 
transition matrix and the input distribution matrix for the random 
forcing function. Several methods have been proposed for the identi- 
fication of the state transition matrix, but no one has proposed a 
method for identifying the input distribution matrix for systems' with 
a single random input. The purpose of this investigation is to propose 
such a method. The problem may be defined as follows: 

Given : 

1. A system whose behavior may be defined by a set of linear, 
constant -coe fficient , differential equations. 

2. A statistical description of the excitation of the system. 

3. A sequence of observations on the state vector. 

4. An estimate of the state transition matrix of the system. 
Prob lem : 

Determine the input distribution matrix of the system. 

Although this investigation is concerned only with systems with 
random forcing functions, the methods described may be extended to 
systems with both random and control inputs. The method of section 
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Ill may be extended by changing the form of the estimated state transi- 
tion matrix as described by Lee [5]; and that of section IV, by including 
the effect of the control input in the filter equations as described by 
Kalman [3]. The problem is restricted to random forcing functions for 
simplicity. 
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II. BASIC MATHEMATICAL MODEL 



The investigation reported here starts with the assumption that 
there exists a set of n linear, constant-coefficient, differential 
equations of the form 

x(t) “ Ax(t) + Bw(t) (2.1) 

which describe the dynamic behavior of a system subjected to a random 
forcing function, w(t). The solution of these equations is 

x(t) - $(t-t^) x(0) + ( § (t -T ) Bw (t )dT (2.2) 

. ^ fc o 

where 

*(t) -X ." 1 Ksl - A)' 1 ]- 

Because the continuous system dynamics will be observed at discrete 
times, it is necessary to describe the system by discrete or difference 
equations. Introducing a sampling device of period T and a zero- 
order hold on the forcing function allows the system equations to be 
written 

Xk + 1 = $(T)x k . + (2.3) 

where Xfc is the system output vector at a time kT, $ is the nxn discrete 
state transition matrix, F(T) is the lxn input distribution matrix and 
w^_ is the sampled and zero-order held input vector at time kT. The 
components of the input vector w^ are assumed to form a gauss ian white 
sequence of zero mean and known variance during the period of system 
operation under investigation. 

Observations made on the system output are assumed to be scalar 
linear combinations of the states of the system which are contaminated 
by additive gaussian white noise of zero mean and known variance. At 
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(2.4) 




where H is the lxn known measurement matrix and v, is the scalar 

k 

additive measurement noise. The system may then be represented by 
the block diagram of Figure 1. 

In the general problem of identification of linear dynamic systems, 
it is necessary to estimate the elements of the state transition matrix 
before proceeding to identify the input distribution matrix. For the 
purposes of this investigation, this estimation is assumed done by 
using the method described by Lee [5]. The estimation of the state 
transition matrix by this method makes it necessary to change the 
system model to the so-called canonic form. Assuming that the system 
is observable and n-identif iable , the difference equation for the 
system model then becomes 



= ** \ + r * ^ 



(2.5) 




where 



Z-k+l " C2S k+l 



C 



H 

H$ . 




( 2 . 6 ) 



$* = C$C 



0 

-1 0 



I 
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r* = cr = b = 

H* = [ 1 0 0 • • • 0 ] 

The matrix is now the estimated state transition matrix obtained by 
using the method of Lee. The canonic system model can then be 
represented by the block diagram of Figure 2. 
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FIGURE 1. BLOCK DIAGRAM OF THE SYSTEM 



v 



k 





FIGURE 2. BLOCK DIAGRAM OF THE 

SYSTEM IN CANONIC FORM 
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III. DERIVATION OF THE IDENTIFICATION ALGORITHM 



Consider a system with canonic state equations of the form 

*k+l " ** + r * \ 

(3.1) 

Z k+1 = H * ^k+1 + V k 

where w^ is assumed to be a scalar random forcing function. The 
z-trans formed state equation becomes 

z£(z) - z_y.(°) = z( z ) + r * w ( z ). (3.2) 

Assuming that the initial conditions ^(0) are zero or that the 
identification process does not start until the initial conditions 
no longer affect the system output, the system may be represented by 
the signal-flow graph of Figure 3. From this signal-flow graph, the 
transfer function Z{z)/w(z) may be written as 




z 2 n-1 n 



z z z 

Since the elements of the matrix are assumed to be known, the 

problem can now be stated as finding the location of the zeros of 

the transfer function and thus the values of b 13 b~ , . b . 

1 2 n 

From the transfer function, the difference equation relating the 

measured output z to the input w may be written as 
K K 
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W(z) 



N 
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FIGURE 3. SIGNAL-FLOW GRAPH OF THE SYSTEM 



(3.4) 



2 k = \ + tp l 2 k-n + cp 2 2 k-^l + -" + < * > „ z k- b l“k-l- (b 2-'< > n b l ) "k-2 

■W2\l b l>Vj 

-•••• (b n-r C,> n b „-2----- <? 3 b l > “k-^l 

- (b -cp b . .-Cp 0 b )w 

n n n-1 2 1 k-n • 

Shifting indices by +1 and using vector notation, Equation (3.4) 



becomes 



J k+1 






d 4- v. 



(3.5) 



where 



1 



'k-n+1 

2 k-n+2 



i J 



eg = 



% 



CPi 



W. 



W. 



W, 



k- 1 
k-n+l 



d = 



(3.6) 
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Because the elements of both the vectors n and d are unknown, a single 
noise element is defined as 

n k -a t T d + v k 0.7) 

where is a scalar and 

Q = d. w. + d 0 w, - + ...+ d w. . - + v. 

K lk 2 k-1 n k-n+1 k 

s d, w, n + d 0 w, + . . . 4* d w, . 0 + v. 



^k+1 U 1 w k+l ' u 2 "k ' ^n vv k-n+2 ' v k+l 



Q = d w + d 0 w + ... + d w. + v. 

it+n 1 k+n 2 k-n+1 n k+1 k+n 



(3.8) 



From Equation (3.8) it is apparent that the noise element is 
correlated with the elements °k+2- •••' “k+n-b ho ” ever ' «k 

and are independent with respect to each other. This property 

of the noise elements is the basis for the identification algorithm. 
As shown in Appendix A, a set of autocorrelation functions for the 
noise elements may be defined such that 



R(0) = E = O w 2 (d^ 2 + d 2 2 + ... + d^ 2 ) + a v 2 



(3.9) 



n-J 



= 1, 2 , . . . , n-1 



k(j) - e - ° w ^ d, d k+j j - i 

R(j) =0 j > n 

where a is the standard deviation of the gaussian white sequence of 
w 

noise inputs w, and a is the standard deviation of the measurement 
r k v 



noise . 



By rearranging Equation (3.5) such that 

fk - ■ 






J k+1 



(3.10) 
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the noise element is now a linear function of the system outputs 



Z k-n+l’ Z k-n+2’ 



• • • y 



z and is therefore known. Hence, the auto- 
kt i. 



correlation functions of equation (3.9) may be found by taking the 



of n simultaneous nonlinear algebraic equations which may be solved 

for the values of d. . d„. d . The values of b,. b«. .... b may 

12 n 12 n 

i 

then be found from equation (3.6). 

The general procedure for the identification of the input distri- 
bution matrix in a linear system with a random forcing function may be 
outlined as follows: 

1. Measure a scalar linear combination of the states of the system 
plus additive noise. 

2. Identify the canonic form of the state transition matrix of 
the system. 

3. Calculate the values of the noise element Q at each sampling 



4. Calculate the autocorrelation functions R(j) for j =0, 1, . .., 

n- 1 . 

5. Solve the simultaneous equations (3.9) for the elements of d. 

6. Solve equation (3.6) for the elements of b. 

Although the above algorithm does provide a means for identifying 
the input distribution matrix, two problems arise in the identification 
process. The first of these involves deciding when to start the 

identification process. Because the autocorrelation functions are 

* 2 

computed by taking the average of the values of > ^\^k+l’ •••» 



until enough time has elapsed so that the sample mean and variance of 



statistics E 




Equation (3.9) now represents a set 



time . 




at a number of sampling instants, this process must not begin 
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the noise input and the measurement noise approach 0, G and 0, G 

w v 

respectively. If the identification process is begun too soon, the 
computed values of _b will be in error from the actual values. 

The other problem in the identification process can be best shown 
by an example. Assuming n = 2 , the autocorrelation functions then 
become 

2 2 2 2 
R (0) “ g (d n + d 0 ) + G 
w 1 2 v 



R(l) - O w 2 d 1 d 2 



(3.11) 



Solution of Equation (3.11) gives 




2,2 

0 

v 



2.2 
a ] 

v 




1/2 




(3.12) 

1/2 



and from Equation (3.6) 



b 



1 



- d 



1 



(3.13) 

• b 2 ~ d 2 + ^2 b l 

It is apparent from Equation (3.12) that the identification of the 
elements b^, b^ involves a decision as to which combination of signs 
to use as the correct one. For higher order systems this decision 
becomes more complex since the number of possible sign combinations 
increases . 
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IV. KALMAN FILTER APPROACH 



An alternative to the identification algorithm of the preceding 
section may be derived from the equations of the Kalman filter 
described below. Consider again a system with the canonic state 
equations of the form 

*k +1 ■ ** + r * 

(4.1) 

Z k+1 ‘ H * *k+l + Vl 



where w^ is a scalar random forcing function. The canonic state 
vector may be estimated by the technique described by Kalman [3] 
where, given k previous measurements, the estimate at the 



sampling instant is given by 



4/k ' ** 4-i/k-i + °k u k • H * ** 4-i/k-i 1 



(4.2) 



th . 



where is the filter weighting or gain applied at the k iteration 
and H* is the lxn known measurement matrix. The gain is calculated 



from a set of recursive equations of the form 

c k - Vk-i h * t ' h * p k/k-i h * t + R r ' 



p k/k = <p - G k H *> p k/k-i 



(4.3) 



Vl/k * ** P k/k ** + o 



where P^^ and are t * ie est iniation and prediction covariances. 



respectively, and R and Q are given by 

A 2 

R = E [v k Z ] 

A 2 T 

q = r* e [w ] r* 1 

K 



(4.4) 
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Since all elements necessary to predict the canonic state vector 
^ are known except the input distribution matrix F*, the system model 
based on this estimation scheme will produce suitable estimates for ^ 
only if the correct T* is used in the filter equations. In order to 
determine if the correct T* is being used, it is necessary to define a 
performance index to compare the performance of the model with that of 
the system. Because the actual system output is available only in the 
form of the measurement , the performance index must be defined in 
terms of and the predicted measurement z^^ ^ as given by 

<4 - 5) 

which is obtained from the system model. 

Before defining the performance index to be used, it is necessary 
to define some other terms. The filter is said to have reached the 
steady-state condition when the prediction covariance matrix does not 
change from one sampling instant to the next. Under the steady-state 
condition the actual prediction measurement variance is given by 

1 K+M a ? 

EZ (M) = - £ (z - z , L ) (4.6) 

i“K 

where K is the value of k at the sampling instant when the filter 
reaches the steady-state condition. The integer M is chosen such 
that the relation 

|EZ(M+1) - EZ(M) I < 0.05|eZ(M)| (4.7) 

is valid. The predicted measurement variance is defined as 

pz - E t< z k -\/k-i> 2) <4 - 8) 
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( 4 . 9 ) 



which, as shown in Appendix B, may be written as 



PZ = H* P H* + R 
K 



where P is the prediction covariance under the steady-state conditions 
of the filter. 

From the above definitions the performance index J is defined as 
the difference between the actual and the predicted measurement 
variances and may be written as 

J = EZ(M) - PZ 



J = 



r - 1 K+M - 

- E (z. - z. . . ) 

M i i/i-l . 



H* H* + R 
K 



( 4 . 10 ) 



Since it is desired that the output of the model ^ an< ^ t ^ e out P ut 

of the system z. be equal for perfect operation of the model, J is to 
K. 

be minimized. Because both EZ(M) and PZ are functions of the input 
distribution matrix T*, the minimum for J will result when the correct 
T* is used in the model. The problem may now be stated for this 
approach as identifying T* so as to minimize the performance index J. 

Because it would waste time and effort, a tria 1-and-error search 
for the correct T* is undesirable. Therefore, it is necessary to 
derive some method of changing the elements of T# used in the model 
after each calculation of J, based on the value of J. The method 
investigated is based on a gradient search where the elements of F* 
are changed by some amount according to the values of a sensitivity 
function and J. The sensitivity function must reflect the change in 
J caused by changing T*. The method of changing T# may be written as 

r* ^ r* + f (j, r*> j (4. n> 

new old — ’ 7 
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where the subscripts "old" and "new" represent the value of T* just 
used in the model and the value to be tried next, respectively. The 
sensitivity function is denoted by _f (J, T*) where f is a vector 
function . 

The problem is now to find a suitable sensitivity function to use. 

Two functions and their failures during simulation are presented in 

the next section. Because the performance index is a function of the 

measurements z , z, r , _ , ..., z Trlx , (which are functions of all past 

K K+l K+M ^ 

measurements) and the predicted measurements during the interval K, 
K+M, (which are also functions of all past measurements), the problem 
of finding a suitable sensitivity function is a complicated one. At 
present no suitable function has been found. 
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V. COMPUTER SIMULATIONS 



In order to determine the accuracy of the proposed identification 
algorithm, a second-order system was simulated on the IBM 360/67 com- 
puter. For the system chosen the canonic state transition matrix, the 
input distribution matrix, and the measurement matrix were 






0 1 

-1/16 1/2 



r* = 



i 

1/2 



H* = jh 0 J 



(5.1) 



In the simulation, the actual and the canonic-state transition matrices 
were assumed equal so that no error in the identification process would 
result from a difference in the two, and also for simplicity. The random 
forcing function and the measurement noise were obtained from a random 
number generator which gives a gaussian white sequence of any desired 
mean and standard deviation. The means and standard deviations for the 
two sequences were 



E [w ] ~ 0 

K. 

E [w k 2 ] « 2 
E [v R ] =0 



(5.2) 



The system dynamics were simulated and the measurements z were used 
to obtain the autocorrelation functions R(0) and R(l)« These values 
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were ‘then used to compute the values of and at each sampling 
instant from the equations 



1/ 2—i -,1/2 



b L (k) = [o. 125 [r(0) - 1 + [ (R(0) - l) 2 - 4R(1) 2 ] ] ] 



(5.3) 



1 / 2 -. -. 1/2 



b 2 (k) = [o . 125 [r( 0) - 1 - [ (R(0) - l) 2 - 4R(1) 2 ] ] ] 



The system and the identification process were simulated for 1000 
sampling intervals of period one second. An ensemble of 200 simu- 
lations was run and an ensemble average for R(0) , R(l), b^(k), and 
b 2 (k) was taken. The results at various sampling intervals are shown 
in Table I. The ensemble averages of the differences between the 
actual values of T* and the computed values are also shown. 

As can be seen from Table I, there is a bias in the autocorrelation 
function R(l) which produces a bias in the estimates for b^ and b^ . 

No explanation for the bias can be given. Simulation of another, 
similar second-order system also produced a bias, equal to -1.0. 

For the algorithm to produce the desired results, it is necessary to 
eliminate the bias. 

A second-order system was also simulated for the Ka lman-f liter 
approach to the identification problem. The values used for the 
simulation were 






0 1 

1/8 1/4 

— 



r* = 



i 

1/2 

— — 



H* = 



L 1 




(5.4) 
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The random forcing function and the measurement noise were generated as 



described above with 



E [wj = 0 
E [w k 2 ] = 3 
E [v fc ] = 0 



(5.5) 



E [v 2 ] = R = 4 . 

K 



The performance index J was calculated for various values of F* with 
the results given in Table IX. It can be seen that J is a minimum 
when the correct T* is used. 

Two sensitivity functions were also calculated during the simulation. 
The first as derived in Appendix C is given by 



As shown in Table II, neither of these sensitivity functions is 
suitable since neither produces a consistent pattern by which to 
change F*. The failure of these functions appears to lie in the fact 
that neither measures the sensitivity of the measurements to the values 

of r*. 




(5.6) 



The second, also derived in Appendix C, is given by 




* (5.7) 
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k 

2 

10 

100 

200 

300 

400 

500 

600 

700 

800 

900 

999 

R(O) 

R(l) 



b^k) 


1-b^k) 


b 2 (k) 


l/2-b 2 (k) 


E ^\ 


e t^k^k+d 


0.0 


0.0 


0.0 


0.0 


5.11 


0.0 


0.84 


0.16 


0.65 


-0.15 


5.25 


-0.31 


1.003 


-0.003 


0.645 


-0.145 


5.17 


-0.49 


1.009 


-0.009 


0.635 


0. 135 


5.18 


-0.50 


1.010 


-0.01 


0.633 


-0.133 


5.17 


-0.50 


1.012 


-0.013 


0.633 


-0.133 


5.18 


-0.50 


1.016 


-0.016 


0.631 


-0. 131 


5.20 


-0.50 


1.017 


-0.017 


0.634 


-0.134 


5.21 


-0.51 


1.016 


-0.016 


0.633 


-0. 133 


5.20 


-0.51 


1.018 


-0.018 


0.634 


-0.135 


5.22 


-0.51 


1.019 


-0.019 


0.636 


-0.136 


5.23 


-0.51 


1.019 


-0.019 


0.635 


-0.135 


5.23 


-0.51 


r 2 (d L 2 + 


. 2, . 2 
d„ ) + a s 
2 v 


2 

“ 5.0 (should match EfO, 

K. 


]) 




2 d l d 2 = 


0.0 (should 


match E 









TABLE I. Results of Simulation 
of Identification Algorithm 
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V 


b 2 


EZ (M) 


PZ 


J 


(J. 


r*> 


i 2 » 


, r*) 


1/2, 


1/4 


14.7 


6.7 


8.7 


-8.4, 


-2.3 


-2.4, 


-0.6 


1, 


1/4 


14.6 


13.4 


1.2 


-17.7, 


-0.3 


12.1, 


0.2 


3/2, 


1/4 


15.0 


24.8 


-9.8 


-26 . 6 , 


0.4 


24.1, 


-0.3 


1/2, 


1/2 


14.1 


7.4 


6.7 


-6.3, 


-4.2 


0.1, 


0.04 


*1, 


1/2 


13.7 


13.3 


0.4 


-18.0, 


1.1 


12.5, 


-0.8 


3/2, 


1/2 


14.3 


24.5 


-11.2 


- 27.3, 


2.1 


24.7, 


-1.9 


1/2, 


3/4 


15.7 


8.5 


7.2 


-0.3, 


-7.6 


0.1, 


1.3 


1, 


3/4 


19.9 


11.0 


8.9 


-60.5, 


60 . 6 


23.0, 


-23.0 


3/2, 


3/4 


15.0 


23.0 


-8.0 


-31.6, 


8.9 


28.1, 


-7.9 



^Correct T* 



TABLE II. Results of Simulation of 
Kalman Filter Approach 
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VI. CONCLUSIONS 



The identification of the input distribution matrix of a linear 
system operating in a stochastic environment is a complicated process, 
but it can be done to some degree by the proposed algorithm. Further 
investigation is needed to determine the reason for the bias in the 
autocorrelation function R(l). A general procedure for choosing the 
proper signs to use in the solution of the nonlinear algebraic equations 
is also desirable. 

The Kalman-f ilter approach would be ai excellent method to use in 
the identification process if a suitable sensitivity function can be 
found. This function must reflect the sensitivity of the performance 
index to the changes in the input distribution matrix. A possible 
first step in finding the function is to examine the sensitivity of 
the measurements to the elements of the matrix. 
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APPENDIX A 



DERIVATION OF THE AUTOCORRELATION FUNCTIONS 
FOR THE NOISE ELEMENT 
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APPENDIX B 



DERIVATION OF THE PREDICTED MEASUREMENT VARIANCE 
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APPENDIX C 



DERIVATIONS OF THE SENSITIVITY FUNCTIONS 

c)PZ 

The first sensitivity function to be derived is during steady- 

state operation of the filter. This is done by finding an expression 
for the prediction covariance during steady-state operation, equating 
elements on both sides of the equation, and then taking the desired 
partial derivatives of the implicit functions. The derivation follows. 
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Another sensitivity function of the form 
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is to be considered now. It may be found by expressing the filter gain 
in steady state as a function of the steady-state prediction covariance. 
The required partial derivatives may then be found as above. The follow- 
ing derivation is done assuming steady-state operation of the filter. 
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where is as given above. 



40 



BIBLIOGRAPHY 



1. Demetry, J. S., Class Notes, EE4414, Naval Postgraduate School, 

Academic year 1968-1969. 

2. Ho, Y. C. and Whalen, B. , "An Approach to the Identification and 

Control of Linear Dynamic Systems with Unknown Parameters," 

IEEE Transactions on Automatic Control , No. 3, July, 1963. 

3. Kalman, R. E., "A New Approach to Linear Filtering and Prediction 

Theory," Transactions of the ASME , V. 82, March, 1960, pp. 

35-45. 

4. Kaplan, W. , Advanced Calculus , pp. 90-97, Addis on-Wes ley Publishing 

Co. , Inc . , 1957 . 

5. Lee, R. C. K. , Optimal Estimation, Identification, and Control , 

pp. 79-126, MIT Press, 1964. 



41 



INITIAL DISTRIBUTION LIST 



No. Copies 

1. Defense Documentation Center 20 

Cameron Station 

Alexandria, Virginia 22314 

2. Library, Code 0212 2 

Naval Postgraduate School 

Monterey, California 93940 

3. Naval Ordnance Systems Command Headquarters 1 

Department of the Navy 

Washington, D. C. 20360 

4. Professor James S. Demetry, Code 52De 1 

Department of Electrical Engineering 

Naval Postgraduate School 
Monterey, California 93940 

5. LTJG Michael R. Polk, USN 1 

1021 Woodlark 

Tyler, Texas 75701 



42 



Security Classification 



DOCUMENT CONTROL DATA - R & D 

(Security classilication 0 / title, body ol abstract and indexing annotation must be entered when the overall reoort Is classilied) 


\ originating ACTIVITY (Corporate author) j 

Naval Postgraduate School 
Monterey, California 93940 


!?0. REPORT SECURITY CLASSIFICATION 

Unc lassi f ied 


2b. GROUP 


3 REPORT TITLE 

Identification of the Input Distribution Matrix for Linear Dynamic Systems 
Operating in a Stochastic Environment 


4 DESCRIPTIVE NOTES (Type of report and.inctusive dates) 

Master's Thesis; June, 1969 


5 AUTHORtSl (First name, middle initial, la st name) 

Michael Raymond Polk 




fi. REPORT DATE 

June, 1969 


7a. TOTAL NO OF PAGES 7b. NO OF REFS 

41 5 


6«. CONTRACT OR GRANT NO. 
6. PROJEC T NO. 

C. 

d. 


9a. ORIGINATOR’S REPORT NUMBERIS) 


9b. OTHER REPORT NO(S) o(/ief numbers that may be assigned 

this report) 


10 DISTRIBUTION STATEMENT 

Distribution of this document is uni 


imited . 



11. SUPPLEMENT ARY NOTES 



12. SPONSORING MILITARY ACTIVITY 



13. ABSTRACT 



An algorithm for the identification of the input distribution matrix 
of a linear, stationary system operating in a stochastic environment is 
derived. The identification is accomplished by defining a set of auto- 
correlation functions for a noise element composed of a linear combination 
of the input distribution matrix elements and the random excitations of 
the system. Another possible identification method employing a Kalman 
filter is discussed and the problems associated with its derivation are 
presented. Results of computer simulations for both methods are included. 



FORM I47O 

1 nov es I “ / w 

S/N 0101 -807-681 1 



(PAGE 1) 



43 



Security Classification 



A-31408 



Security Classification 



KEY WO R OS 



Identification 
Estimation 
Stochastic 
Discrete Systems 



DD ,"'“,,1473 'back, 

S/N 01 01 -807-682 1 





LIN 


K A 


LINK B 


LINK C 




ROLE 


W T 


ROLE 


W T 


ROLE 


W T 

















44 



Security Classification 



A- 3 I 409 



lhesP683 

Identification of the input distnbu on 




3 2768 000 99271 3 

DUDLEY KNOX LIBRARY 



